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ABSTRACT 


Fluid transients, also known as water hammer, can have a significant impact on the 
design and operation of both spacecraft and launch vehicle propulsion systems. These transients 
often occur at system activation and shutdown. The pressure rise due to sudden opening and 
closing of valves of propulsion feed lines can cause serious damage during activation and 
shutdown of propulsion systems. During activation (valve opening) and shutdown (valve closing), 
pressure surges must be predicted accurately to ensure structural integrity of the propulsion 
system fluid network. 

In the current work, a network flow simulation software (Generalized Fluid System 
Simulation Program) based on Finite Volume Method has been used to predict the pressure 
surges in the feed line due to both valve closing and valve opening using two separate 
geometrical configurations. The valve opening pressure surge results are compared with 
experimental data available in the literature and the numerical results compared very well within 
reasonable accuracy (< 5%) for a wide range of inlet-to-initial pressure ratios. A Fast Fourier 
Transform is preformed on the pressure oscillations to predict the various modal frequencies of 
the pressure wave. 

The shutdown problem, i.e. valve closing problem, the simulation results are compared 
with the results of Method of Characteristics. Most rocket engines experience a longitudinal 
acceleration, known as “pogo” during the later stage of engine burn. In the shutdown example 
problem, an accumulator has been used in the feed system to demonstrate the “pogo” mitigation 
effects in the feed system of propellant. The simulation results using GFSSP compared very well 
with the results of Method of Characteristics. 
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INTRODUCTION 


Fluid transients (also known as Water hammer) have significant impact in design and 
operation of spacecraft and launch vehicle propulsion systems. The pressure rise due to sudden 
opening and closing of valves of propulsion feed line can cause serious damage during activation 
and shutdown of propulsion systems. Pressure surge occurs when either a propellant feed line 
system is opened or closed suddenly by using control valves. The accurate prediction of pressure 
surge is quite important from structural integrity point of view of the propulsion systems. Most 
rocket engines experience a longitudinal instability, known as “Pogo” during the later stage of 
engine burn. To mitigate pogo oscillation, an accumulator is used in the feed system of liquid 
propellant to detune the vibration coupling of the feed line with the vehicle structural frequencies. 

There have been numerous studies to predict the pressure surge in pipeline-reservoir 
systems. Hearn [1] has studied the effect of valve opening in propellant loading line from a 
reservoir both analytically and experimentally. The Method Of Characteristic (MOC) is one of the 
most widely used semi-analytical method for water hammer surge prediction [2, 3]. This method 
is based on solving ordinary differential equation along line of characteristics. MOC, however, is 
not particularly suited for typical fluid network with branching, solid-to-fluid heat transfer and 
phase change. Majumdar and Flachbart [4] have used the network flow analysis simulation 
based on finite volume method to compute the fluid transients of long cryogenic pipeline and they 
have compared their results with that of the Method of Characteristic simulation results. 

In most recent study, Lee [5] has studied the fluid transient problem in predicting the 
pressure rise in entrapped air in a horizontal pipe experimentally. This problem is an example of 
sudden opening of valve and sudden acceleration of liquid mass into the air space. The presence 
of air in the pipe causes excessive pressure rise in the pipe. 

In the present study, two separate problems representing pressure rise due to valve 
opening and valve closing are solved numerically using the finite volume based network fluid 
analysis simulation program (Generalized Fluid System Simulation Program) and the results are 
compared with the experimental results of Lee [5] and semi-analytical MOC results respectively. 

A Fast Fourier Transform analysis has been conducted to convert the pressure-time 
characteristics into a frequency domain to predict the frequency band of the pressure oscillations. 

PROBLEM DESCRIPTION 

Two separate problems are considered for the study. The first problem is an example of a 
valve opening problem where liquid is allowed to enter a pipe line having entrapped air by 
opening a ball valve. The second problem is an example of sudden valve closing. These two 
problems are described in brief below. 

a) Valve Opening Problem 

A long pipe is attached to a reservoir at one end containing liquid water and closed at the 
other end as shown in the figure 1 . The liquid water and entrapped air regions in the pipe are 
separated by a ball valve. The dimension of the pipe and other controlling parameters such as 
reservoir-to-air pressure ratio, length of air column etc are taken identical to that used by Lee [5] 
so that the numerical results can be compared to the experimental data of Lee [5]. The ball valve 
is opened from a 0% opening to 100% opening by controlling the angle of the ball valve and this 
is shown in figure 2. The reservoir pressure is a few times higher than the pressure of the 
entrapped air (air is assumed to be at atmospheric pressure). The ratio of reservoir pressure to 
the initial pressure (P R ) varies in the range of 2 to 7, i.e. the reservoir pressure (p R ) range being 
29.4 psi to 102.9 psi. The dimension of the pipe and other parameters are identical to that used 
by Lee [5] in his experimental work, so that the current study results can be compared with Lee’s 
experimental results. Apart from the initial pressure ratio, another controlling parameter is the 
ratio of initial length of the entrapped air column to the total length of the pipe (a = L g /L T ). The 
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initial length for the water volume in the pipe (L|) is fixed to 20 ft, and initial length of air column in 
the pipe (L g ) varies from a low of 1 .23 ft to 1 6.23 ft, the value of a ranging from 0.0579 to 0.448 
respectively. The pipe diameter is 1.025 inches. The entrapped air and water is initially at 14.7 
psia and 60 F. 



The ball valve does not open until about 0.15 sec, and gradually starts opening and it 
opens 100% at about 0.4 second. Figure 2 shows the ball valve angle position with time and 0 
degree refers to full closed and 90 degree refers to full open position. In the present numerical 
model this is accounted for by providing valve area change history. 



Time (Seconds) 

Figure 2. Bail Valve Angle Change with Time [5] 
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b) Sudden Valve Closure 

A long pipe line carrying liquid oxygen is subjected to sudden valve closure. An 
accumulator containing helium gas is placed upstream of the valve to reduce the amplitude and 
frequency of the pressure oscillation. A schematic diagram of the system is as shown in figure 3. 
The pipe is 400 ft long and 1/4 inch diameter. The accumulator is placed just before the valve. 



Accumulator 


j Valve position 


•* 




400 ft 


Figure 3. Schematic diagram of the valve closing problem with accumulator. 

In the present study, the effect of the accumulator is studied by removing it from the feed 
line and placing it back just before the valve. It is assumed that the accumulator is charged with 
gaseous helium, and represents a pogo suppressing device in the feed line. 


RESULTS AND DISCUSSION 

In this section the numerical method based on finite volume principle will be described in 
brief followed by the results for two separate problems. The generalized fluid simulation software 
program (GFSSP) is based on the numerical method described and instrumental in solving both 
the problems of valve opening and closing. The results of the valve opening problem are 
compared with the experimental data of Lee [5] and Lee and Martin [6] and the results of the 
valve closing problem are compared with the Method of Characteristics. 

NUMERICAL MODELING 


GFSSP is a finite volume based network analysis software [7], which resolves a given flow 
network into a finite number of nodes interconnected by branches. The nodes consist of a) 
Boundary Nodes where pressure, temperature and species concentrations are specified, and b) 
Internal Nodes where pressure, temperature and species concentrations are calculated by 
solving mass, energy and species conservation equations respectively. The momentum 
equations are solved in the branches to calculate the flowrate entering or leaving the node. The 
flow could be steady or unsteady, compressible or incompressible, and with or without heat 
transfer. Thermodynamic properties at each node are computed by GASP [8] or WASP [9]. The 
system of equations is solved by a combination of Newton-Raphson and successive substitution 
method. GFSSP has a user friendly Graphical User Interface, VTASC that allows user to develop 
GFSSP model of complex network by point and click paradigm. 


GFSSP MODELS: 
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a) Numerical Model for the Valve Open Problem: 

The physical problem shown in Figure 1 is numerically modeled by dividing the region 
from the exit of the reservoir to the valve by using a total of 12 nodes, out of which node 1 is the 
boundary node as shown in Figure 3. The restriction is used to control the flow of water into the 
pipe. The volume of node 12 is adjusted as the volume of water expands and the corresponding 
surrounding air volume reduces. This is done to adjust mutual change in volume of the water 
column and the entrapped air. 
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Figure 4. Finite Volume Model of the Flow Network 

This is done in GFSSP through the user subroutine by adding a volume source for the working 
fluid (water in this case). This volume source is computed based on the relation given below. 

Change in node volume: 

As soon as the valve is opened, the water in the pipe is going to compress the entrapped 
air, and this volume increase of water is adjusted as a volume source for node 12. The air volume 
would reduce by the same amount. Using ideal gas assumption and thermodynamic relation, the 
volume source for the working medium (water) is given as: 

A V wa t cr AVair — m a i r R an -T 1 2 Ap/p*" 

Where, Ap = change in pressure at node 12 with time = p - p , where both p and p represent 
pressure computed at node 12, at the current and previous time steps respectively. 


Computation of Air pressure: 

The entrapped air is considered stationary and so the momentum and other conservation 
equations are not solved for air. The bulk pressure in the air domain is computed by using 
thermodynamic relations. As the whole process is very fast (valve opening is about 4 seconds), 
the compression in the air is assumed to be adiabatic. 


Pair 
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Additional adjustment for momentum source at node 12: 

The pressure developed within the air will try to prevent water column movement and this 
resistance needs to be incorporated in the model and this is done as a momentum sink to the last 
control volume (branch 11-12) of the computation domain. 

Momentum Source for branch 1112 = -(p air - pi 2 )A 
where, A is the cross-sectional area. 
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b) Numerical Model for the Valve Closure Problem 
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Figure 5. GFSSP Model of Valve Closure Problem 

Figure 5 shows a flow network consisting of a 400 feet pipeline modeled with five internal 
nodes (2, 3, 4, 5 and 6) and five (12, 23, 34, 45 and 56) branches. Branch 67 represents a valve 
which closes in 0.1 second. Node 1 and 7 represents boundary nodes representing supply and 
outlet pressures respectively. Pressures and temperatures are computed in the internal nodes 
by solving the mass and energy conservation equations and flowrate are computed in branches 
by solving the momentum conservation equations. In this simulation the accumulator was 
modeled in node 5 as a node with variable volume. Part of node volume is occupied by Helium 
gas; the gas volume changes due to compressibility. As the pressure wave travels along the 
pipeline, the gas volume fluctuates to provide compliance to the system. 

RESULTS FOR VALVE OPEN CASE 


The computed results of GFSSP is compared with that of Lee and Martin [6] for three 
different geometrical configurations: (a) a = 0.448 (b) a = 0.195 and (c) a = 0.058. The table 
below shows the values of air co lumn length as a function of alph a. 


a 

*-g 

L t 

0.4480 

16.23 

36.23 

0.1952 

4.85 

24.85 

0.0579 

1.23 

21.23 


Table 1: List of Geometrical Parameters for Numerical Case Study 


GRID INDEPENDENCE TEST: A grid independence study was done for the numerical solution to 
see the effect of grid size. For this purpose, a case study was done with the inlet to initial 
pressure ratio P R = 7 and air length to total length ratio a = 0.448. Figure 6(a) shows that the 
entrapped air pressure does not change appreciably when the number of nodes is doubled. 
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Figure 6(a). Grid Independence study of numerical solution on air pressure 
TIMESTEP INDEPENDENCE: 

In order to get a time step independent solution of the problem, the simulation was 
carried out with time steps of 0.01 , 0.005 and 0.0025 sec. These time steps satisfy the Courant 
condition of Courant No. = UAx/Ax where U is characteristic velocity, Ax is the time step and Ax 
is the spatial interval. CFL condition needs that the Courant no. should be less or equal to 1 for 
a stable solution. Figure 6(b) shows that a reduction of time step from a time step of 0.005 s does 
not change the numerical results considerably and hence for the entire simulation, a time step of 
0.005 sec is adequate for reasonably accurate solution. 



Figure 6(b). Time independence study of the numerical solution. 
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VALIDATION WITH EXPERIMENTS: 

The numerical results of GFSSP are compared with the experimental data of Lee [5] for 
comparing the pressure developed in the entrapped air for three different geometrical 
configurations (a) a = 0.445 with Pressure ratio (Pr) of 4 and 7, (b) a = 0.1952 with Pressure 
ratio (P R ) of 4 and 7 and (c) a = 0.058 with Pressure ratio (P R ) of 3 and 5. For same pressure 
ratio, the maximum pressure rise will be for smaller a. Figure 7 show the plot of air pressure 
variation with time for case (a) a = 0.448 with two different pressure ratios (I) P R = 4 and (II) P R = 
7. 
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Figure 7. Validation of air pressure for (I) P R = 4 and (II) P R = 7 at about 45% initial air 

volume (a = 0.45) 
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Figure 8. Air Pressure at various Pressure ratios for case (b), a= 0.2 (20% air) 



Figure 9. Air Pressure at various Pressure ratios for case (c), a= 0.06 (6% air) 

Figures 8 and 9 show the pressure rise in the entrapped air when relative less air is 
present in the entrapped pipe initially (approximately 20% and 6% by volume respectively). The 
air pressure rises to approximately 260 psia and 380 psia respectively when the reservoir 
pressure is 73.5 psia for 20% and 6% initial air respectively. The comparison with the 
experimental data (represented by dotted lines in the figures) show quite good agreement as far 
as amplitude of the pressure oscillations are concerned even at higher reservoir pressure, and 
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the agreement is excellent (within 5% or so) at medium and low reservoir pressure (44 psia or 
less). The discrepancy at higher pressure could be due to the rigid pipe assumption and due to 
not accounting for pipe deformation. 


WATER PRESSURE AND FAST FOURIER TRANSFORM 

The pressures in the water column at different nodes (spatial locations) are plotted as a 
function of time as shown in figure 1 0 below. Node 2 is the close to the entry of pipe from the 
reservoir and node 1 1 is almost at the end of the pipe line. The pressure oscillation is due to the 
water hammer effect, and also due to compressed air pushing the water column back from 
expanding. The frequency of pressure oscillations can be better analyzed by converting the 
transient pressure response to a frequency domain as shown in Figure 1 1 . 
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Figure 10. Pressure Transients as various spatial locations 
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Figure 11. Nodal Pressure response in Frequency Domain. 

Figure 1 1 shows frequency band of the primary, secondary and other higher order 
modes. This transformation is being done by using a First Fourier Transformation of the pressure- 
time response using the following transformation function. 
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X(e“)= ]Tx[nK' t!n 

n=-oo 

Here, x[n] represents the time domain function and X[Q] in the frequency domain. The 
real part of the complex field represents the amplitude and the imaginary part gives the 
frequency. The lower frequency represents the primary mode oscillation frequency and the 
higher ones for the other modes. 


RESULTS FOR VALVE CLOSURE PROBLEM 


Figure 12 shows the GFSSP simulation results for the valve closing problem (no 
accumulator) with two different grids (20 and 40) and results are compared with the Method of 
Characteristic results. Further increase in number of grid points to 80 nodes did not change the 
numerical results appreciably. The amplitude and frequency matched quite well between the 
network flow analysis (GFSSP) simulation with the semi-analytical Method of Characteristic 
results. 
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Figure 12: Comparison of numerical results of GFSSP simulation with MOC. 

Figure 13 shows the effect of having accumulator in the pipe line, and it has been shown 
that the accumulator dampens the pressure pick and the pressure oscillation eventually 
completely zeroed down. 
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Sudden Valve Closing: GFSSP Results with and w/o Accumulator 



Figure 13: Effect of Accumulator on Pressure Rise 


SUMMARY AND CONCLUSIONS 

This paper demonstrates that the finite volume based network flow analysis method 
implemented in GFSSP can accurately predict fluid transient during rapid opening and closing of 
valve. It can also model the compliance caused due to the compressibility of the gas. The 
simulation results from GFSSP compared reasonably well with the experimental data for the 
valve opening problem and with that of method of characteristics for the valve closing problem. 
For the valve opening problem with entrapped air in the pipe, the agreement of computed results 
with the experimental data of Lee [5] is excellent at high ratio of air to water, and in particularly 
when the ratio of inlet pressure to initial air pressure is low (<= 4). The deviation between 
experiment and prediction at higher pressure ratio is due to the assumption of rigid pipe while the 
experiments were carried out in plastic tubes that may have added to the compliance. Future 
work will include fluid-structure interaction as well as the effect of heat transfer and phase change 
in fluid transients. 
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